/*BORWEIN'S ALGORITHM*/

/*
Comando de compilacao:
gcc T1_BorweinP.c -o BorweinP -lgmp -pthread

Grupo 5, turma B, 2013
Programacao Concorrente

Algoritmo:

a = 6 - sqrt(32)
y = sqrt(2) - 1
ynext = (1-(1-y⁴)^(1/4)) / (1+(1-y⁴)^(1/4))
anext = a*(1+ynext)^4 - 2^(2*iteracao+3)*ynext*(1+ynext+(ynext)^2)

anext vai convergir quarticamente em 1/pi, portanto,
o numero de digitos corretos aproximadamente quadruplica
a cada iteracao. Como o numero de digitos corretos apos
a virgula na primeira iteracao é 7, necessitamos de 11
iteracoes para obtermos um numero superior a 10milhoes 
de digitos corretos.
*/

#include <stdlib.h>
#include <stdio.h>
#include <pthread.h>
#include <gmp.h>
#include <time.h>

#define ITER 11 /*numero de iteracoes*/
#define PRECISION 33300000 /*precisao*/

/*Declarando variaveis com alta precisao*/
mpf_t a, anext;
mpf_t y, ynext;
mpf_t pival;
mpf_t auxa1, auxa3, auxy1, auxy2;

/*fator que soluciona problema em representar uma equacao
explicada dentro do while do codigo*/
int fator = 8;

/*funcoes para substituir codigo e colocar dentro de threads*/

void * dothread1 (void * args)
{
	/*a = 6 - sqrt(32)*/
	mpf_set_str (anext, "32", 10);
	mpf_sqrt (anext, anext);
	mpf_ui_sub (anext, 6, anext);
	return NULL;
}

void * dothread2 (void * args)
{
	/*y = sqrt(2) - 1*/
	mpf_set_str (ynext, "2", 10);
	mpf_sqrt (ynext, ynext);
	mpf_sub_ui (ynext, ynext, 1);
}

void * dothread3 (void * args)
{
	/*auxy1*/
	mpf_pow_ui (auxy1, y, 4);
	mpf_ui_sub (auxy1, 1, auxy1);
	mpf_sqrt (auxy1, auxy1);
	mpf_sqrt (auxy1, auxy1);
	mpf_ui_sub (auxy1, 1, auxy1);
	return NULL;
}

void * dothread4 (void * args)
{
	/*auxy2*/
	mpf_pow_ui (auxy2, y, 4);
	mpf_ui_sub (auxy2, 1, auxy2);
	mpf_sqrt (auxy2, auxy2);
	mpf_sqrt (auxy2, auxy2);
	mpf_add_ui (auxy2, auxy2, 1);
	return NULL;
}

void * dothread5 (void * args)
{
	/*auxa3*/
	mpf_pow_ui (auxa3, ynext, 2);
	mpf_add (auxa3, ynext, auxa3);
	mpf_add_ui (auxa3, auxa3, 1);
	mpf_mul (auxa3, ynext, auxa3);
	mpf_mul_ui (auxa3, auxa3, fator);
	fator = fator*4;
	return NULL;
}

void * dothread6 (void * args)
{
	/*auxa1*/
	mpf_add_ui (auxa1, ynext, 1);
	mpf_pow_ui (auxa1, auxa1, 4);
	mpf_mul (auxa1, a, auxa1);
	return NULL;
}



/*funcao principal*/

int main()
{
	/*variaveis para controlar o tempo de execucao*/
	time_t time_start;
	time_t time_end;
	double time_dif;
	
	/*pegar o tempo inicial*/
	time(&time_start);
	
	/*Setando a precisao*/
	/*10000000 digitos representam 10^10000000 numeros distintos. x bits representam 2^x numeros distintos, portanto temos aproximadamente x = 33.219.000 bits*/
	mpf_set_default_prec (PRECISION);
	
	/*variavel para controlar o numero de bits a serem comparados*/
	unsigned long int cmp_bits;
	
	/*variavel de controle de iteracoes*/
	int iteracoes;
	/*variaveis para lidar com arquivo*/
	char* arquivo;
	FILE * F;
	
	
	/*Threads*/
	pthread_t thread1, thread2, thread3, thread4, thread5, thread6;
	
	/*Inicializando variaveis*/
	mpf_init (a);
	mpf_init (anext);
	mpf_init (y);
	mpf_init (ynext);
	mpf_init (pival);
	mpf_init (auxa1);
	mpf_init (auxa3);
	mpf_init (auxy1);
	mpf_init (auxy2);
	iteracoes = 0;
	
	arquivo = "pivalueparallel.txt";
	
	/*Atribuindo valor inicial para "a" e "y"*/
	/*a = 6 - sqrt(32)*/
	pthread_create(&thread1, NULL, dothread1, NULL);
	/*y = sqrt(2) - 1*/
	pthread_create(&thread2, NULL, dothread2, NULL);
	/*espera thread1 e thread2 acabarem*/
	pthread_join(thread1, NULL);
	pthread_join(thread2, NULL);
	
	
	/*Calcula iteracoes*/
	while (iteracoes < ITER)
	{
		/*Guarda as variaveis antigas para serem usadas*/
		mpf_swap (a, anext);
		mpf_swap (y, ynext);
		/*ynext = (1-(1-y⁴)^(1/4)) / (1+(1-y⁴)^(1/4))*/
		/*auxy1*/
		pthread_create(&thread3, NULL, dothread3, NULL);
		/*auxy2*/
		pthread_create(&thread4, NULL, dothread4, NULL);
		/*espera thread3 e thread4 acabarem*/
		pthread_join(thread3, NULL);
		pthread_join(thread4, NULL);
		/*ynext = (auxy1) / (auxy2)*/
		mpf_div (ynext, auxy1, auxy2);
		
		/*anext = a*(1+ynext)^4 - 2^(2*iteracao+3)*ynext*(1+ynext+(ynext)^2)*/
		/*anext = auxa1-auxa2*auxa3*/
		/*auxa3*/
		pthread_create(&thread5, NULL, dothread5, NULL);
		/*auxa1*/
		pthread_create(&thread6, NULL, dothread6, NULL);
		/*espera thread5 e thread6 acabarem*/
		pthread_join(thread5, NULL);
		pthread_join(thread6, NULL);
		/*anext = auxa1-auxa3*/
		mpf_sub (anext, auxa1, auxa3);
		
		/*pi = 1 / anext*/
		mpf_ui_div (pival, 1, anext);
		
		/*imprime na tela o valor de pi com algumas casas*/
		gmp_printf("pi = %.30Ff\n", pival);
		
		/*Incrementa contador de iteracoes*/
		iteracoes++;
	}
	
	/*pegar o tempo de termino*/
	time(&time_end);
	
	/*fazendo a diferenca para pegar o tempo corrido*/
	time_dif = difftime(time_end, time_start);
	
	/*Salvando o valor do pi em 10milhoes de casas em arquivo*/
	F = fopen (arquivo, "w");
	mpf_out_str(F,10,10000010, pival);
	fclose(F);
	
	/*liberando espaço utilizado*/
	mpf_clear(pival);
	mpf_clear(a);
	mpf_clear(anext);
	mpf_clear(y);
	mpf_clear(ynext);
	mpf_clear(auxa1);
	mpf_clear(auxa3);
	mpf_clear(auxy1);
	mpf_clear(auxy2);
	
	
	
	/*imprime tempo corrido*/
	printf("Tempo de execução = %.0f segundos\n", time_dif);
	
	return 0;
}
